Notebook 5: Rational Method and Travel Times¶
In this notebook, we look at two ways of estimating a runoff hydrograph from precipitation data.
First, we’ll use the rational method to approximate peak flow and the maximum water level at the basin outlet, and then we’ll use an open-source library to make a higher resolution estimate of flow accumulation paths and stream networks using a digital elevation model.
# import required packages
import pandas as pd
import numpy as np
import math
# advanced statistics library
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.patches as mp
import matplotlib.colors as colors
# SEE COMMENTS ABOUT PYSHEDS LIBRARY IN NEXT CELL
from pysheds.grid import Grid
import warnings
warnings.filterwarnings('ignore')
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
from bokeh.models import LinearColorMapper, LogTicker, ColorBar, BasicTickFormatter
from bokeh.io import output_notebook
output_notebook()
%matplotlib inline
Import Precipitation Data¶
For this exercise, we will use historical climate data from the Meteorological Service of Canada (MSC) station at Whistler, BC.
# calibration data
df = pd.read_csv('../data/Whistler_348_climate.csv',
index_col='Date/Time', parse_dates=True)
# note that the 'head' command shows the first five rows of data,
# but in this case the columns are abbreviated.
# print(df.head())
# list all the columns
# print('')
# print('__________')
for c in df.columns:
print(c)
stn_name = df['Station Name'].values[0]
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
<ipython-input-2-7fd4599e7107> in <module>
1 # calibration data
----> 2 df = pd.read_csv('../data/Whistler_348_climate.csv',
3 index_col='Date/Time', parse_dates=True)
4 # note that the 'head' command shows the first five rows of data,
5 # but in this case the columns are abbreviated.
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options)
608 kwds.update(kwds_defaults)
609
--> 610 return _read(filepath_or_buffer, kwds)
611
612
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in _read(filepath_or_buffer, kwds)
460
461 # Create the parser.
--> 462 parser = TextFileReader(filepath_or_buffer, **kwds)
463
464 if chunksize or iterator:
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in __init__(self, f, engine, **kwds)
817 self.options["has_index_names"] = kwds["has_index_names"]
818
--> 819 self._engine = self._make_engine(self.engine)
820
821 def close(self):
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in _make_engine(self, engine)
1048 )
1049 # error: Too many arguments for "ParserBase"
-> 1050 return mapping[engine](self.f, **self.options) # type: ignore[call-arg]
1051
1052 def _failover_to_python(self):
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in __init__(self, src, **kwds)
1865
1866 # open handles
-> 1867 self._open_handles(src, kwds)
1868 assert self.handles is not None
1869 for key in ("storage_options", "encoding", "memory_map", "compression"):
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/parsers.py in _open_handles(self, src, kwds)
1360 Let the readers open IOHanldes after they are done with their potential raises.
1361 """
-> 1362 self.handles = get_handle(
1363 src,
1364 "r",
~/Documents/code/data_analysis/lib/python3.8/site-packages/pandas/io/common.py in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
640 errors = "replace"
641 # Encoding
--> 642 handle = open(
643 handle,
644 ioargs.mode,
FileNotFoundError: [Errno 2] No such file or directory: '../data/Whistler_348_climate.csv'
Plot the measurement¶
It’s a good idea to visualize the data we’re working with.
# plot flow at Stave vs. precip at the closest climate stations
p = figure(width=900, height=400, x_axis_type='datetime')
p.line(df.index, df['Total Precip (mm)'], alpha=0.8,
legend_label="Total Precip [mm]", color='dodgerblue')
p.line(df.index, df['Total Snow (cm)'], alpha=0.8,
legend_label="Total Snow [cm]", color="firebrick")
p.line(df.index, df['Total Rain (mm)'], alpha=0.8,
legend_label="Total Rain [mm]", color='green')
p.legend.location = 'top_left'
p.legend.click_policy = 'hide'
p.xaxis.axis_label = 'Date'#
p.yaxis.axis_label = 'Daily Rain [mm] / Snow[cm] Volume'
show(p)
Simplified Version of Rainfall-Runoff¶
First, isolate a single precipitation event to use for estimating a runoff hydrograph. Let’s find a nice week for skiing:
fig, ax = plt.subplots(1, 1, figsize=(16,4))
sample_start = pd.to_datetime('2014-12-01')
sample_end = pd.to_datetime('2014-12-15')
sample_df = df[(df.index > sample_start) & (df.index < sample_end)][['Total Precip (mm)', 'Total Snow (cm)', 'Total Rain (mm)']]
# print(sample_df.head())
ax.plot(sample_df.index, sample_df['Total Precip (mm)'], label="Total Precip [mm]")
ax.plot(sample_df.index, sample_df['Total Snow (cm)'], label="Total Snow [cm]")
ax.plot(sample_df.index, sample_df['Total Rain (mm)'], label="Total Rain [mm]")
ax.set_xlabel('Date')
ax.set_ylabel('Precipitation')
ax.set_title('{}'.format(stn_name))
plt.legend()
<matplotlib.legend.Legend at 0x7fefd8b33c40>
First, imagine we are some unfortunate parking lot attendant working a shift in Whistler Village at Parking Lot 5, and we are told by our cruel supervisor we have stand at the lowest point of the parking lot: a catchment basin with an area of \(1 km^2\) where water runs off into FitzSimmons Creek. The sky looks angry, but we’re running late for work and put on our running shoes instead of our sturdy waterproof boots.
Next, assume the travel time is effectively zero across our entire basin (precipitation takes no time to travel to the outlet once it falls on the parking lot surface). Is this a reasonable assumption in general?
Under these assumptions, lets reconstruct a runoff hydrograph at the outlet. First, look at the precipitation data over the twelve days of the big storm.
print(sample_df)
Total Precip (mm) Total Snow (cm) Total Rain (mm)
Date/Time
2014-12-02 0.0 0.0 0.0
2014-12-03 0.0 0.0 0.0
2014-12-04 4.4 6.0 0.0
2014-12-05 12.0 0.0 6.0
2014-12-06 12.0 0.0 6.0
2014-12-07 2.8 0.0 2.8
2014-12-08 30.5 0.0 30.5
2014-12-09 89.6 0.0 89.6
2014-12-10 74.4 0.0 74.4
2014-12-11 19.2 0.0 19.2
2014-12-12 3.6 0.0 1.8
2014-12-13 0.0 0.0 0.0
2014-12-14 0.0 0.0 0.0
Convert Volume to volmeteric flow units¶
Runoff is typically measured in \(\frac{m^3}{s}\), so convert \(\frac{mm}{day}\) precipitation to \(\frac{m^3}{s}\) runoff.
# convert to runoff volume
drainage_area = 1 # km^2
# runoff is typically measured in m^3/s (cms for short -- cubic metres per second),
# so express the runoff in cms
sample_df['runoff_cms'] = sample_df['Total Rain (mm)'] / 86.4
print(sample_df)
Total Precip (mm) Total Snow (cm) Total Rain (mm) runoff_cms
Date/Time
2014-12-02 0.0 0.0 0.0 0.000000
2014-12-03 0.0 0.0 0.0 0.000000
2014-12-04 4.4 6.0 0.0 0.000000
2014-12-05 12.0 0.0 6.0 0.069444
2014-12-06 12.0 0.0 6.0 0.069444
2014-12-07 2.8 0.0 2.8 0.032407
2014-12-08 30.5 0.0 30.5 0.353009
2014-12-09 89.6 0.0 89.6 1.037037
2014-12-10 74.4 0.0 74.4 0.861111
2014-12-11 19.2 0.0 19.2 0.222222
2014-12-12 3.6 0.0 1.8 0.020833
2014-12-13 0.0 0.0 0.0 0.000000
2014-12-14 0.0 0.0 0.0 0.000000
If the channel outlet has a rectangular shape of width 2m, how tall should our boots be? Assume a 2% slope, and find a reasonable assumption for the roughness of asphalt.
Recall the Manning equation:
Where:
n is the manning roughness
A is cross sectional area of the flow
R hydraulic radius (area / wetted perimeter)
S is the channel slope
w_channel = 1.5 # m
S = 0.005 # channel slope
n_factor = 0.017 # rough asphalt
def calc_Q(d, w, S, n):
"""
Calculate flow from the Manning equation.
"""
A = d * w
wp = w + 2 * d # wetted perimeter
R = A / wp
return (1/n) * A * R**(2/3) * S**(1/2)
def solve_depth(w, n_factor, Q, S):
"""
Given a flow, a roughness factor, a channel slope, and a channel width,
calculate flow depth.
"""
e = 1 / 100 # solve within 1%
d = 0
Q_est = 0
n = 0
while (abs(Q_est - Q) > e) & (n < 1000):
Q_est = calc_Q(d, w, S, n_factor)
# print(Q, Q_est, abs(Q_est - Q))
d += 0.001
n += 1
# print('solved in {} iterations'.format(n))
return d
# For each timestep, we want to solve for the depth of water at our outlet
sample_df['flow_depth_m'] = sample_df['runoff_cms'].apply(lambda x: solve_depth(w_channel, n_factor, x, S))
plt.plot(sample_df.index, sample_df['flow_depth_m'])
plt.ylabel('Flow depth [m]')
Text(0, 0.5, 'Flow depth [m]')
Not only are our feet wet, but if we happen to be there the peak it’s potentially dangerous. As little as 10-15cm of water moving fast enough can sweep you off your feet.¶

More Complex Implementation¶
As discussed in class, precipitation takes time to travel from where it fell to the basin outlet. Next we will estimate the runoff response in a real catchment, just upstream from the parking lot example in the FitzSimmons Creek basin.
Step 1: Instantiate a grid from a DEM raster¶
Some sample data is already included, but for extra data, see the USGS hydrosheds project.
# grid = Grid.from_raster('data/n45w125_con_grid/n45w125_con/n45w125_con', data_name='dem')
grid = Grid.from_ascii(path='../data/notebook_5_data/n49w1235_con_grid.asc',
data_name='dem')
# reset the nodata from -32768 so it doesn't throw off the
# DEM plot
grid.nodata = 0
# store the extents of the map
map_extents = grid.extent
min_x, max_x, min_y, max_y = map_extents
Plot the DEM¶
NOTE: The cell below may take up to 30 seconds to load. Please be patient, it is thinking really hard.
The code below will plot the Digital Elevation Model (DEM).
Do you recognize any features of the terrain? Can you locate where it is?
Hover over the map (or touch if using a touchscreen) to see the coordinates in decimal degree units.
What does the precision of the coordinates represent?
i.e. what does 5 decimal places in decimal degrees equate to in kilometers?
You can interact with the plot by using the tools on the left (in vertical order from top to bottom):
pan: move around the map
box zoom: draw a square to zoom in on
wheel zoom: use the mousewheel (or pinch gesture on a touchscreen) to zoom in
box zoom: draw a square to zoom in on
tap: not yet implemented (but you can see the coordinates)
refresh: reset the map
hover: see the coordinates when hovering over the map with a mouse or pointer
# set bokeh plot tools
tools = "pan,wheel_zoom,box_zoom,reset,tap"
# show the precision of the decimal coordinates
# in the plot to 5 decimal places
TOOLTIPS = [
("(x,y)", "($x{1.11111}, $y{1.11111})"),
]
# create a figure, setting the x and y ranges to the appropriate data bounds
p1 = figure(title="DEM of the Lower Mainland of BC. Hover to get coordintes.", plot_width=600, plot_height=int(400),
x_range = map_extents[:2], y_range = map_extents[2:],
tools=tools, tooltips=TOOLTIPS)
# map elevation to a colour spectrum
color_mapper = LinearColorMapper(palette="Magma256", low=-200, high=2400)
# np.flipud flips the image data on a vertical axis
adjusted_img = [np.flipud(grid.dem)]
p1.image(image=adjusted_img,
x=[min_x], # lower left x coord
y=[min_y], # lower left y coord
dw=[max_x-min_x], # *data space* width of image
dh=[max_y-min_y], # *data space* height of image
color_mapper=color_mapper
)
color_bar = ColorBar(color_mapper=color_mapper, #ticker=Ticker(),
label_standoff=12, border_line_color=None, location=(0,0))
p1.add_layout(color_bar, 'right')
show(p1)
# -123.15512, 49.41293
#-123.14657, 49.41080
-123.14350, 49.40251